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Abstract 

We investigate the behavior of impurity fields immersed in a larger condensate field in 1, 2, and 3 di- 
mensions. We discuss the localization of a single impurity field within a condensate and note the effects of 
surface energy. We derive the functional form of the attractive interaction between two impurities due to 
mediation from the condensate. Generalizing the analysis to N impurity fields, we show that within various 
parameter regimes a crystal of impurity fields can form spontaneously in the condensate. Finally, the system 
of condensate and crystallized impurity structure is shown to have nonclassical rotational inertia, which is 
characteristic of superfluidity, i.e. the system can be seen to exhibit supersolid behavior. 
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Introduction - Recent observations in solid helium [1] have generated enormous speculation 
about the existence of supersolidity, first predicted almost forty years ago [2]. This issue remains 
controversial, partly due to the difficulty of investigating the helium system both theoretically 
and experimentally. Trapped atomic Bose-Einstein condensates (BECs), which exhibit superfluid 
behavior, are a promising alternative system to probe supersolidity as one can not only tune the 
atomic interactions, but also calculate many of their properties from first principles. In this paper 
we present a model of small impurity fields immersed within a larger condensate (impurities in 
Bose liquids have been considered in many contexts [3]) that mimics some of the fundamental 
properties of supersolids, namely that it breaks the translational symmetry of the Hamiltonian 
as well as exhibits off-diagonal long range order as evidenced by nonclassical rotational inertia 
(NCRI). Unlike previous supersolid proposals in dilute BECs I4D, here the crystal scale emerges 
spontaneously from the system rather than being externally imposed. 

The present model is realizable with current experimental technology. For instance, the dis- 
tinct impurity fields can be produced by utilizing isotopes or by transferring atoms from a larger 
condensate via Raman pulses into different atomic levels in an atomic trap. One can also explore 
parameter space by tuning the coupling constants via a magnetic and optical Feshbach resonance. 

The model- In this paper, we consider a large BEC denoted by the field ip coupled to TV small 
distinguishable impurity fields denoted by Xk- To permit an uncluttered description of the system's 
nontrivial properties, we shall assume that all impurities interact with the same coupling constants. 
(We briefly consider more realistic systems in the discussion.) Our system is therefore governed 
by the following Hamiltonian: 

r l~i i N N 1 N N 

H= / ^|V^ + iH 4 + AH 2 ^| Xfc | 2 + ^^-|Vx fc | 2 + |El^| 4 + 7E^| 2 |^l 

'' |_ k=l k=l k=l j<k 

where A is the coupling of the condensate to the impurity fields, 70 is the self-interaction of the im- 
purity fields, and 7 is the interaction between impurity fields. If 70 = this Hamiltonian describes 
distinguishable impurities that do not self-interact. We shall restrict ourselves to positive param- 
eters. In dilute BECs these coupling constants are directly proportional to the atomic scattering 
length with a proportionality constant of 4:7rh 2 /m. As mentioned above, these scattering lengths 
are in principle tunable via magnetic fields in current experiments. 

From eq.®, the system's dynamics are governed by TV + 1 coupled nonlinear Schrodinger 



da(l) 
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equations: 

N 



id t *lj = + \^ + \^\Xk\ 2 (2) 

k=i 



idtXk = -^Axk + 7o\Xk\ 2 Xk + 7Xk^T\Xj\ 2 + M^\ 2 Xk (3) 

where A is the Laplacian in D spatial dimensions. In this system there is particle conservation of 
the large condensate field N = J \i[j\ 2 d D x and of each impurity field n k = J \xk\ 2 d D x, and we 
assume n k <C N. The total energy (OQ) and the total linear momentum P = Im J ifr* d D x + 
I m ^2k=i m I Xt^Xk d D x are also conserved. 

Instabilities of the uniform state and collapse of the system - Let us consider the nontrivial 
structures that impurity fields can generate in a condensate. To determine the critical point at 
which these structures emerge, we begin with the uniform miscible state, which is stable for certain 
values of the coupling constants. Ignoring surface tension effects (discussed below), one can show 
that the uniform state is stable if and only if the (TV + 1) x (N + 1) matrix, M, defined by 
the quadratic form of the potential energy in (OQ), is positive semidefinite, i.e. if all eigenvalues 
of M are positive [5]. One can also prove that if M is negative semidefinite the system will 
experience finite-time collapse [6]. Furthermore, because the density is non-negative, one can 
extend the criteria in [6] and show that the system would also experience a finite-time collapse if 
M is conegative or positive subdefinite. We will focus mainly on the regime where the uniform 
state has a modulational instability (implying that at least one eigenvalue of M is negative) and, 
since the coupling constants are assumed positive, we do not consider situations proven to exhibit 
finite-time collapse. 

The N + 1 eigenvalues of M are 70 — 7, which is TV — 1 degenerate, and |(1 + (N — 1)7 + 



70) ± \\/{N — 1)7 + 70 — l) 2 + 47VA 2 . Therefore, a modulational instability occurs if 7 > 70 

or if A > - — -j= , highlighting the difference between the system's two distinct types of 

modulational instabilities. In the first regime, 7 > 70, each impurity localizes individually whether 
or not the system is phase separated from the condensate. In the second regime, the condensate 
phase separates from the impurity fields. 

A single impurity: Interpolating between self-localization and phase separation - We now 
discuss the system of a small impurity field embedded in and interacting with a large condensate. 
A variational argument will show that there is a critical value of the coupling parameter (dependent 
on surface tension effects) between the condensate and impurity beyond which the impurity self- 



localizes. 

The ground-state wave functions in 1, 2, and 3 dimensions should possess the following char- 
acteristics: that x(r) be localized such that x(r) —> as r —> oc, and its mass be fixed such that 

D_ 

nk = C D J |x/c| 2 ^ D_1 ^, where C D = |^ry is the surface of a unit sphere in D spatial dimensions; 
that there be a depletion where the impurity is positioned, and ip(r) —> ip = cte as r — > oc. The 
energy to be minimized is the Hamiltonian (OQ) for one impurity relative to the energy of uniform 
state assuming i/j(r) Re- 
consider the real, normalized trial functions of the form x(r) = \^^^^^f{ ar ) anc * 
R(r) = Ro (1 — a x( r ) 2 )' where a and a > are variational parameters determined by mini- 
mizing the energy estimates (note that any real finite system does not allow a —> 0) and f(r) is 
localized such that f(r) — > as r —> oc. The normalization constant is Af D = Jq 00 f{x) 2 x D ~ 1 dx. 
With a formal expansion for the energy, one can easily show that a ^ in the limit a —> . 
The energy in D dimensions is thus bounded by 

E < E{a) = e a 2 + e x a D + e 2 a 2+D + e 3 a 3D (4) 

where the constants e , e 2 , and e 3 are all positive numbers (detailed calculations will be presented 
elsewhere), but ei = fci (70 — A 2 ) n\ (where fci is positive)may change its sign if A 2 > 70. 

It is important to emphasize that the variational approach only gives us an upper bound on the 
ground-state energy. If for a nonzero a the lowest variational energy is negative (recall the energy 
for the nearly uniform state is positive and approaches zero as a —> 0), then we know the uniform 
state is unstable and can be reasonably sure that self-localization has occurred. However, if the 
lowest variational energy is non-negative, we cannot determine whether or not self-localization 
occurs. Nonetheless, this expression of energy as a function of the variational parameter a, eq. 
©, does give useful insight into the localized impurity solution in 1-, 2-, or 3-dimensional space. 
When 61 > 0, the energy is a monotonically increasing function of a, implying that the nearly 
uniform ground state (a —> 0) minimizes the energy. 

For D = 1, the dominant term at small a is 61 a. Thus when 61 is negative, i.e. A 2 > 70, a 
supercritical transition occurs from a homogeneous state to a localized impurity state. 

For D = 2, the first and second terms are of the same order. As is the case for D = 1, 
there is a second-order transition towards a localized impurity state if A 2 > 70 + with k\ = 
4C2 A/2 r °oc -ft vl , . The instability of the homogeneous state is shifted from the bulk condition 

Jo J\ x ) x ® x 

described above. This shift has a simple interpretation: the presence of a 1 /m factor means that 
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this term comes from the kinetic energy of the impurity J ^ | Vxk | 2 > so the shift is created by the 

f°° fix) 2 xdx 

curvature a of the localized structure -it is a kind of surface tension. Furthermore, jog f( x ^4 x ' ^ 
the ratio of an interface energy to a bulk energy, bears the hallmark of a surface-tension effect. 

For D = 3, the situation is more subtle. As discussed above, the variational argument states 
that a localized solution exists if the minimum energy in © is negative for some critical value, 
a c . A necessary (but not sufficient) condition for a negative-energy ground state is A 2 > 70. We 
can make a more precise estimate as follows: because the energy turns negative for some range 
of a > and the energy expansion grows quadratically near a « (since e > in ©), it is 
sufficient (though not necessary) that E(a c ) is a minimum and that E(a c ) < 0. The critical line 
defining onset of self-localization is given by E(a c ) = and E\a c ) = 0; on one side we can say 
that a localized structure exists, on the other side we cannot be certain. This line maybe written in 
a parametric way by the following: 6/3 7 + 2/3 3 = e^je^ & 7/? 6 + 3/? 2 + £ 1 ^/s 3 2 /2 = 0, 
where (3 = a c (ss/s2) 1 ^ 6 is the parametrization. In the gaussian approximation one gets the upper 
bound sup I (^) 3/4 ^=^, m 3/5 ( ^ o) 2/ 5 } in the case of 7o = 0, which is about 35% higher 

than the numerical result [|7[]. 

Condensate-mediated attraction- Here we discuss the interaction between TV dilute localized 
impurity fields and derive the effect of the perturbed background condensate on the interaction. 
We will show that the interaction has an attractive tail mediated by the condensate in addition to 
the hard-core repulsion arising from the repulsion between impurities. 

Let us consider TV self-localized impurity fields (having satisfied the conditions given in the 
previous section) that weakly modify the uniform condensate, i.e. ip = ip + ipi(r). This as- 
sumption is valid if the impurity is either only weakly localized or sufficiently distant. The latter 
situation is of particular relevance here as we are interested in deriving the long-range attractive 
tail of the interaction. In this approximation, the condensate wave function is linear and may be 
solved with the aid of a Green's function in D-spatial dimensions. The solution is well known and 
can be written explicitly as 

11 e -2^ \x-x f \ inD = 1 
K (2xp \x -x'\) inD = 2 (5) 
-2xIjq\x-x'\ . 
T \ — inD = 3 
\x—x' I 

where K is the modified Bessel's function. We can thus determine ipi(x). Assuming the con- 
densate to be only weakly modified, eliminating the terms that do not depend explicitly on the 
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impurity fields, and explicitly omitting the self-interaction energy of the impurities which is given 
by 170, the interacting energy of the system of TV impurities becomes 

N 



i^k 

where 



E =\Y, j U{x-x')\ Xi {x')\ 2 \xk{x)\ 2 d D xd D x' 



(6) 



4\2 /2 

U(x t - x k ) = ^ D) { Xi - x k ) - —P^G^ixi - x k ). (7) 

Let us consider localized impurities (a\x k — Xi\ 3> 1). One may approximate the k-th impu- 
rity field by \xk(&)\ 2 = n k 8^ D \x — x k ) where x k is the position of the k-th impurity. In this 
case E = \ Y^i^k nin k U(\xi — x k \), where the interaction potential between the two impurities 
is given by Uix) = U(x) of ©. (The yukawa attractive tale for similar systems in 3D were 
discussed in [8]) The first 5-interacting energy term in © is only a crude estimation. The full 
expression © says that a bound state exists but the equilibrium distance is zero. In reality the 
equilibrium distance is on the order of the size of the impurity. Introducing the trial function 
Xfc(x) = V^\/^^ J ? (a\x — x k \) with f(s) = e" s2 , into the interaction energy ©, for D = 3 
on finds U(\x\) = 7^372 e~ a2 ^ 2 — ^°- e 2 ^ lxl . The second term was arrived at by assuming 
a ^> 2^o, i.e. that the impurity is more localized than the Yukawa interaction range 1 / (2^ ). The 
equilibrium distances come from this energy. 

N impurities — crystallization - We now discuss the system in which we have TV interacting 
impurity fields within a much larger condensate. We showed in the previous section that within 
a certain parameter regime the localized impurities exhibit a tunable hard-core attractive weak 
interaction in any D with like particles, a system that has been studied in the context of many 
diverse physical systems (for 3D classical particles see e.g. [9]). We will show that these impurities 
crystallize in two regimes of condensate-impurity interaction. FigureQ]is a phase diagram showing 
four distinct regimes. In phase I, the condensate and impurities are miscible. In phase //, the 
condensate phase separates from the impurities, which remain together and miscible with each 
other in a bubble. No crystallization occurs in either phase. Crystallization does occur in phases 
/// and TV. In phase the impurities remain immersed in the condensate; in phase IV they phase 
separate, forming a crystal within a bubble within the condensate. 

In phase the impurities are attracted through the interaction energy ©. However, there is a 
repulsive hard core, and an equilibrium distance is expected. The system crystallizes as is clearly 
shown through numerics (see Fig. [2]). 
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FIG. 1: Phase diagram. The vertical line represents 7 = 70 while the curve / represents v . All 

numerical simulations of Eqns. (12131) are in two space dimensions in a 64 x 64 periodic plane. The initial 
conditions were uniform miscible states plus small fluctuations. In each sector the left image plots the 
condensate density \ip\ 2 while the right, plots J2k=i \Xk\ 2 - In HI the initial masses are / \^\ 2 d 2 x = 4096 
and J \xk\ 2 d 2 & = 40.96, A = 0.2, 7=1 and 70 = while in IV the initial masses are J \ip\ 2 d 2 x = 4096 
and J \xk\ 2 d 2 x = 40.96, A = 2, 7 = 1.5 and 70 = 0. Note that, as discussed in the text, the surface tension 
effects will shift the boundaries. 



In phase TV, there is phase separation of the condensate from the impurities such that the im- 
purities are confined in a kind of bubble within which the condensate density is near zero. At the 
same time, inside the bubble there is also phase separation of the impurity fields from each other. 
The interaction between the impurities is no longer has an attractive tail because ip « in (|5]). 
The condensate essentially plays the role of container since any impurity that begins to stray from 
the ensemble becomes more strongly attracted to it as more condensate seeps in between, caus- 
ing an attractive interaction between the stray impurity and the ensemble. Finally, crystallization 
arises because of the "pressure" of the condensate boundary and the short-distance repulsion, see 
Fig. WV. 

In some sense, the system in phase IV is at a special limit when the condensate is not cou- 
pled to the impurity fields. Therefore, the case A = deserves special attention. The condensate 
can be described by the pure nonlinear Schrodinger equation, which is quite well understood. 
The impurity fields, however, evolve. Impurity fields that overlap store potential energy and will 
tend to repel each other; all the more so as 7 increases. If the impurity fields \xk\ 2 are local- 
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FIG. 2: Numerical simulation of Eqns. ()2I3|) for N = 36 impurity fields in a 128 x 128 periodic plane. 
Here the initial masses are / \^\ 2 d 2 x = 16384 and / \xk\ 2 d 2 x = 200, A = 0.4, 7=1 and 70 = 0. Left a), 
plots the condensate density \ip\ 2 while the right, b) plots J2k=i \Xk\ 2 - Inhomogeneous regions, i.e. where 
miscible states coexist with the crystalline state, indicate that the system has not yet reached an equilibrium 
state. 

ized about R k with a size 5 then, in the limit 7 ^ 70, the potential energy is dominated by 
2 ^2i^k I \Xi( x )\ 2 \Xk(x)\ 2 d D x. Minimizing the overlap will minimize the potential energy. In 
some sense, the energy may be approximated by the superposition of two-body repulsive inter- 
actions. In a finite box there is an upper bound on the separation between the impurities. One 
may invoke a close-packing argument to find the crystal structure that minimizes the energy of the 
impurity ensemble, sustained only by the external pressure from the condensate at the boundaries. 

The periodic case where all TV impurity fields have the same mass n* provides an interesting, 
solvable example. One can use the minimization approach to determine S as a function of the large 
parameter 7. To begin, let us assume that the impurity field vanishes exactly outside the ball of 
radius 5, that there is no overlap, and that only the nearest neighbor(s) affect the interaction energy. 
The minimum of the total energy H = Ylk=i I 2k \^Xk \ 2 dx can be determined from the solutions 
of the Hemholtz equation —^Axk = £Xk for a ball Vs(R k ) with a Dirichlet boundary condition 
on dVs(Rk) [10]. One can perturb this to determine the effect of a slight overlap. The variational 
parameter 5 will determine the optimal configuration that balances the kinetic energy of a pulse 
and the interaction energy due to a small overlap of impurity fields. The one-dimensional case in 
a periodic domain is a solvable example that we shall present elsewhere. 

Nonclassical rotational inertia- We have shown that in a certain parameter regime impurities 
immersed in a condensate crystallize. We will now show that the crystallized phases-/// & IV still 
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behaves as a superfluid (i.e. it has NCRI), and thus exhibits supersolid characteristics. 

Following 111 111 , consider a cylindrical system of volume V, rotating uniformly about its primary 
axis of rotation e,\. For a small angular rotation the system will possess an energy S = S + AS 
where S is the ground state energy and AS = \l eff uj 2 , I eff being the effective and measurable 
moment of inertia tensor around the e\ axis. Deviations of this tensor from the rigid body rotation 
tensor I RB are called the nonclassical rotational inertia fraction. 

Under rotation, the phases of the condensate and impurity fields are no longer uniform and 
the increase of energy is AS = \ J (^Po(x) |V0| 2 + J2k=i mP&WI^^I 2 ) dx where i/j(x) = 
y/ po(x)e^, Xk(x) = a/ pk(x)e l( ^ k are the nonuniform density of the respective ground states, and 
the phases {(f), <p k } satisfy the no matter flux boundary conditions: h • V {0, (j)k\ = • (ei x r). 
That AS is a quadratic form in uu is evident since (f) and (j) k must be proportional to uu due to this 
boundary condition. The explicit form of I efi and the pre-factor needs deeper consideration (see 
Hi for details). The minimization of AS leads to continuity equations and boundary conditions 
that may be solved using the method called homogenization, which splits cleanly the large (system 
size) and small (impurity size) scales and provides effective average quantities. However, it is not 
possible to obtain a closed expression of in terms of the local density po(x) and Pk(x). 
Nevetherless, it is possible to show that < I RB and that I eff is proportional to the superfluid 
density g ss . In one dimension, the problem is exactly solved — a direct calculation leads to 
Leggett's formula fl] : ^ = (i f v J^) + m £f =1 (^ v ^) _1 . 

Because the impurity fields are localized, their wave functions decay fast in space so the contri- 
bution to the superfluid density from the impurity fields becomes negligible. Therefore, to leading 
order the supersolid properties are due to the large condensate; while modulation by the impurities 
causes the NCRI fraction to be less than unity, the impurities themselves do not contibute to the 
NCRI. 

Discussion - Although we have considered systems with identical impurity interaction con- 
stants, it is interesting to consider a more realistic system of many distinguishable impurities (per- 
haps different atomic levels in a single atom) with different coupling constants. Assuming the 
system is stable and possesses a modulation instability (with the constraints on coupling constants 
discussed above), we expect this system to break the translational symmetry of the Hamiltonian 
as systems discussed in this paper and to possess a NCRI. However, rather than forming a regular 
periodic crystal, we expect it to form a giant coherent molecule or amorphous solid with properties 
directly related to the interaction constants. 
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